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ABSTRACT 

Context. Plasma accreting onto classical T Tauri stars (CTTS) is believed to impact the stellar surface at free-fall velocities, generating 
a shock. Current time-dependent models describing accretion shocks in CTTSs are one-dimensional, assuming that the plasma moves 
and transport energy only along magnetic field lines (yS <k 1). 

Aims. We investigate the stability and dynamics of accretion shocks in CTTSs, considering the case of /? > 1 in the post-shock region. 
In these cases the ID approximation is not valid and a multi-dimensional MHD approach is necessary. 

Methods. We model an accretion stream propagating through the atmosphere of a CTTS and impacting onto its chromosphere, by 
performing 2D axisymmetric MHD simulations. The model takes into account the stellar magnetic field, the gravity, the radiative 
cooling, and the thermal conduction (including the effects of heat flux saturation). 

Results. The dynamics and stability of the accretion shock strongly depends on the plasma/?. In the case of shocks with /J > 10, violent 
outflows of shock-heated material (and possibly MHD waves) are generated at the base of the accretion column and strongly perturb 
the surrounding stellar atmosphere and the accretion column itself (modifying, therefore, the dynamics of the shock). In shocks with 
p ~ I, the post-shock region is efficiently confined by the magnetic field. The shock oscillations induced by cooling instability are 
strongly influenced by /?: for /3 > 10, the oscillations may be rapidly dumped by the magnetic field, approaching a quasi-stationary 
state, or may be chaotic with no obvious periodicity due to perturbation of the stream induced by the post-shock plasma itself; for 
yS ~ 1 the oscillations are quasi-periodic, although their amplitude is smaller and the frequency higher than those predicted by ID 
models. 
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1. Introduction 
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In the last few years, high resolution {R — 600) X-ray observa- 
tions of several young stars accreting material from their cir- 
cumstellar disk (TW Hya, BP Tau, V4046 Sgr, MP Mus and 
RU Lupi) suggested the presence of X-ray emission at tempera- 
ture T - 2 - 5 MK from plasma denser than hh = 10" 



cm 



(iKastner et al, 


2002 


; ISchmitt et al. 2005: Gunther et alj l2006t 


Arsiroffi et al. 


2007 


; Ikobrade & Schmitt ,2007£ ,Argiroffi et al.1 



2009h . The emitting pla sma is too dense to be enclosed inside 



coronal loop structures (Te sta et aLl l2004'). Several authors sug- 
gested, therefore, that this soft X-ray emission component could 
be produced by the material accreting onto the star surface, flow- 
ing along the magnetic field lines of the nearly dipolar stellar 
magnetosphere, and heated to temperature s of few MK by a 
shock at the base o f the accretion column (ICalvet & GuUbriiigl 
[19981 inmzin|[T998h . 

The idea that an accreting flow impacting onto a stellar sur- 
face leads to a shocked slab of material heated at millions de- 
grees is not new. For instance, in the context of degenerate stars, 
several authors have studied the dynamics and energetic of ac- 
cretion shocks and the efifects of radiation on the formation, 
structure, and stability of the shocks (Langeret al. 1981^ 



These studies have shown that the role of thermal and dynami- 
cal instabilities is critical to our ability to model radiative shock 
waves. 

The first attempt to interpret the evidence of soft X-ray emis- 
sion from dense plasma in Classical T Tauri stars (CTTSs) in 
terms of accretion shocks has quantitatively demonstrated that 
some non coro nal features of the X-ray observations of TW Hya 
(iGiintheretal .1.200 7) and MP Mus (Argiroffietal 2007) could 
be explained through a simplified steady-flow shock model. 
However, steady models are known to be not a good approx- 
imation for radiative shocks, since they neglect the important 
effects of local thermal instabilities as well as global shock 



oscillatio ns induced by radiative cooling ([Sutherland & Dopita 
1993; D opita & Suflierlandll 19961; ISafieill 19981; [Sutherland et al 



Chevalier & Imamura 1982; lmamura il985[ iChanmugam et al.l 
1985). Huieirat & Papaloizou ( 1998) have also investigated the 
role of the magnetic field in confining the post-shock accreting 
material and in determining the evolution of the accretion shock. 



I2003ai.k iMignonell2005l) . 

In fact, the first time-dependent ID models of radiative ac- 
cretion shocks in CTTSs have shown quasi -periodic oscillations 
of the sho ck position induced by cooling ( Koldoba et al.ll2008t 
I Sacco et al. 2008). In particular, Sacco et al. (2008) have devel- 
oped a ID hydrodynamic model including, among other impor- 
tant physical eff'ects, a well tested detailed description of the stel- 
lar chromosphere and investigated the role of the chromosphere 
in determining the position and the thickness of the shocked re- 
gion. For h ydrodynamic al simulations based on the parameters 
of MP Mus, ISacco et aTl (l2008h synthesized the high resolution 
X-ray spectrum, as it would be observed with the Reflection 
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Grating Spectrometers (RGS) on board the XMM-Newton satel- 
lite, and found an excellent agreement with the observations. 

Up to date time-dependent models of accretion shocks in 
CTTSs have been one dimensional, assuming that the plasma 
moves and transports energy only along magnetic field lines. 
This hypothesis is justified if the plasma yS •« 1 (where fi - gas 
pressure / magnetic pressure) in the shock-heated material. The 
photospheric magne tic field magnitude of C TTSs is believed to 
be around 1 kG (e.g. lJohns-KruU et al.ll999l) . Such a strong stel- 
lar field is enough to confine efliciently accretion shocks with 
particle density below 10'"^ cm"^ and temperature around 5 MK 
{p < 0.3). However, recent polarimetric measurements indicate 
that, in some cases, the photospheric magnetic field strength 
could be less than 200 G dValenti & Johns-Krull|[2004l) and the 
plasma yS in the slab may be around 1 or even larger. In these 
cases, the magnetic field configuration in the post-shock region 
may change, influencing the physical structure of the material 
emitting in X-rays and the stability of the accretion shock. 

The low-y6 approx imation was challenged by recent findings 
of iDrake et al.l ( l2009l) . In case of y6 <«c 1, the radiative shock in- 
stability is expected to lead to detectable periodic modulation of 
the X-ray emission from the shock-heated plasma, if the den- 
sity and velocity of the accretion stream do not change over 
the time interval considered (in agreement with predictions of 
ID models). However, the analysis of soft X-ray observations 
of the CTTS TW Hya (whose emission is believed to arise pre- 
dominantl y from accretion sh ocks) produced no evident periodic 
variations dPrake et al.l 12009 ). These authors concluded, there- 
fore, that ID models may be too simple to describe the multi- 
dimensional structure of the shock and that the magnetic field 
may play an important role through the generation and damping 
of MHD waves. 

In this paper, we investigate the stability and dynamics of ac- 
cretion shocks in cases for which the low-y8 approximation can- 
not be applied. We analyze the role of the stellar magnetic field 
in the dynamics and confinement of the slab of shock-heated 
material. To this end, we model an accretion stream propagating 
through the magnetized atmosphere of a CTTS and impacting 
onto its chromosphere, using 2D axisymmetric MHD simula- 
tions and, therefore, an explicit description of the ambient mag- 
netic field. We investigate cases of yS > 1 in the post-shock region 
for which the deviations from ID models are expected to be the 
largest. In Sect. |2] we describe the MHD model, and the numeri- 
cal setup; in Sect.[3]we describe the results; in Sect.|4]we discuss 
the implications of our results and, finally, in Sect. |5] we draw 
our conclusions. 



2. MHD modeling 

The model describes an accretion stream impacting onto the sur- 
face of a CTTS. We assume that the accretion occurs along mag- 
netic field lines linking the circumstellar disk to the star and con- 
sider only the portion of the stellar atmosphere close to the star. 

The fluid is assumed to be fully ionized with a ratio of spe- 
cific heats y - 5/3. The model takes into account the stellar 
magnetic field, the gravity, the radiative cooling, and the ther- 
mal conduction (including the effects of heat flux saturation). 
Since the magnetic Reynolds number » 1 considering the typi- 
cal velocity (10^ cm s"') and length scale (lO'' cm) of the system, 
the flow is treated as an ideal MHD plasma. The impact of the 
accretion stream is modeled by solving numerically the time- 




Temperature [Kj 



Fig. 1. Radiative losses for an op tically thin plasrri a from the 
APED VI. 3 atomic line database jSmith et alj|200lh, assuming 



the metal abundances 0.5 the solar values jAnders & Grevessa 
09891). 



dependent MHD equations (written in non-dimensional conser- 
vative form): 



dp 

dt 



+ V ■{pu)^Q, 



dpu 

— + V ■ (puu - BB) + VP. = pg . 

ot 



^ + V ■ [uipE + P,) - B(u ■ B)] = 

pu-g-V -F^- n^riiiAiT) , 

dB 

— + V ■ (mB - Bm) = , 

ot 



(1) 
(2) 



(3) 



(4) 



where 



P, = P + 



e + -\? -H 



1 B2 

2 p 



are the total pressure, and the total gas energy (internal energy, 
e, kinetic energy, and magnetic energy) respectively, t is the 
time, p - finiHrtii is the mass density, p - 1.28 is the mean 
atomic mass (assurning m etal abundances 0.5 the solar values; 
lAnders & Grevessd 19891) . niu is the mass of the hydrogen atom, 
«H is the hydrogen number density, u is the gas velocity, g is the 
gravity, T is the temperature, B is the magnetic field, is the 
conductive flux, and A(T) represents the optically thin radia- 
tive losses per unit emission measur43 (see Fig. [TJ derived with 
the PINTof ALE spectral code (Kashyjip_& Drake 2000) with the 
APED VI. 3 atomic line database dSmith et al.ii2001) . assuming 
the same metal abunda nces as before (as de duced from X-ray 
observations of CTTSs; iTelleschi et al.ll2007l) . We use the ideal 
gas law, P - (j - l)pe. 

The thermal conductivity in an organized magnetic field is 
known to be highly anisotropic and it can be highly reduced in 
the direction transverse to the field. The thermal flux, therefore. 



' Note that the radiative losses are dominated by emission hnes in 
the temperature regime common to CTTSs and are set equal to zero for 
T < 10'* K. 
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Fig. 2. Initial geometry of the system in cylindrical coordinates. 
The stellar surface lies on the r axis and the unperturbed stellar 
magnetic field is uniform and oriented along the z axis (vertical 
lines). The accretion stream propagates downwards through the 
stellar corona with velocity Mstr- The z-axis is the symmetry axis 
of the problem. 



is locally split into two components, along and across the mag- 
netic field lines, Fc = F\\ i + F± j, where (see, for instance, 
lOrlando et all2008l) 



Fu = 



1 1 
+ 



\[qspi]\\ [q'satlii 



(5) 



F^^\ 



1 1 

+ 



spijj. 



[qs. 



to allow for a smooth transition between the classical and satu- 
rated conduction regime. In Eqs.|5j [^spiln ™d [^spi]± represent 
the classical conductive flux along and across the magnetic field 
lines (Hiitzer 1962) 



[^spilii = -^ii[vr]|| ^ -9.2 X lo-'T'/^ [vr] 



[qspi]± = -/fx[Vr]x ^ -3.3 X 10-"^;;;^ [WT]^ 



(6) 



7^1/252 



where [Vr]|| and [Vr]j^ are the thermal gradients along and 
across the magnetic field, and k\\ and (in units of erg s"' 
K"' cm"') are the thermal conduction coefficients along and 
across the magnetic field lines, respectively. The saturated flux 
along and across the ma gnetic field lines, [^satln and [^satlj., are 
(ICowie & McKedlI977l) 



[q'satlii = -sign([Vr]||) 54>pcl, 



(7) 



[g'satjj 



-sign([Vr]^) 5cf>pct 



where Cs is the isothermal sound speed, and (p is a number of 
the order of unity; we se t (p = I according to the values sug- 
gested for stellar coronae ( Giulianill984l:lBorkowski et aljl989l; 
[padevev et al. 2002, and references therein). 

We solve the MHD equations using cylindrical coordinates 
in the plane (r, z), assuming axisymmetry and the stellar surface 
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Fig. 3. Initial hydrogen number density (dashed line) and tem- 
perature (dotted line) as a function of height above the stellar 
surface, z, for the unperturbed stellar atmosphere. 

lying on the r axis (see Fig. |2|. The initial unperturbed stel- 
lar atmosphere is assumed to be magneto-static and to con- 
sist of a hot (maximum temperature ^ 10^ K) and tenuous 
(hh ~ 2 X 10^ cm"^) corona linked through a steep transition 
region to an isothermal chromospher^ that is uniformly at tem- 
perature 10"* K and is 8.5 x 10'^ cm thick (see Fig. [3]). The choice 
of an isothermal chromosphere is adopted for ease of model- 
ing, and diff'erent choices of more detailed chromospheric mod- 
els have been shown not to lead to significantly different results 
(Sacco et al. 2009, in preparation). The unperturbed stellar mag- 
netic field is uniform, oriented along the z axis and perpendicu- 
lar to the stellar surface. The gravity is calculated assuming the 
star mass M - 1 .2Mq and the star radiu s - 1.3j^(T) whic h 
is appropriate for the CTTS MP Mus (see Argi roffi et al]|2007l) . 
Different choices of stellar parameters should not lead to differ- 
ent results. 

The accretion stream is assumed to be constant and propa- 
gates along the z axis. Initially the stream is in pressure equilib- 
rium with the stellar corona and has a circular cross-section with 
radius r^a = 5x10'^ cm; its radial density distribution is given 
by 



n,str('-) 



cosh [o- (r/r,tX] 



(8) 



where nstrO and Wcoi are the hydrogen number density at the 
stream center and in the corona, and cr = 20 is a steepness 
parameter The above distribution describes a thin shear layer 
(~ 0.1 Tstr) around the stream that smooths the stream density to 
the value of the surrounding medium and maintains numerical 
stability, allowing for a smooth transition of the Alfven speed at 
the interface between the stream and the stellar atmosphere. The 
stream temperature is determined by the pressure balance across 
the stream boundary. 

The velocity of the stream along the z axis also has a radial 
profile: 



MiitrO 



cosh [o- (r/rstr)"^] 



(9) 



where MstrO is the free fall velocity at the center of the stream. 
The smooth shear layer of velocity avoids the growth of random 



- Note that the radiative losses are set to zero in the chromosphere to 
keep it in equilibrium. 
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Table 1. Relevant parameters of the simulations 



Name of 


Ifiol [G] 


Geometry 


Mesh 


time 


simulation 








covered [s] 


HD-ID 




Cartes. ID 


1024 


3000 


By-01 


1 


cylind. 2D 


512 X 1024 


3000 


By- 10 


10 


cylind. 2D 


512 X 1024 


4000 


By-50 


50 


cylind. 2D 


512 X 1024 


3000 


By-01 -HR 


1 


cylind. 2D 


1024 X 2048 


1000 


By- 10-HR 


10 


cylind. 2D 


1024 X 2048 


1000 


By-50-HR 


50 


cylind. 2D 


1024 X 2048 


2000 



perturbations with wavelengths of the order of the grid size (e.g. 
[Bodo et al. 1994). 

The simulations presented here describe the evolution of the 
system for a time interval of about 3000 s. We used the accretion 
parameters (veloc ity and density) that match the soft X-ray emis- 
sion of MP Mus ( Argiroffi et al.l2007l) . namely nstro = 10" cm"^ 
and MstrO - -500 km s~' at a height z = 2.4 x 10'° cm above the 
stellar surfac^ The initial field strengths of \Bq\ - 1, 10, 50 G 
in the unperturbed stellar atmosphere (see Table [T]l maintain 
> 1 in the post-shock region, where yS = P/iB^/^n) is the ra- 
tio of thermal to magnetic pressure. For instance, \Bo\ = 10 G 
corresponds to /3 ranging between » lO'*, at the base of the chro- 
mosphere, and ^ 1Q~^ up in corona. We also performed a ID hy- 
drodynamic simulation in order to provide a baseline for the 2D 
calculations. This simulation corresponds to the strong magnetic 
field limit (yS ^ 1), in which the plasma moves and trans ports en- 
ergy e xclusively along the magnetic field lines (see iSacco et al.l 
120081) . 

The calculatio ns described i n this paper were performed us- 
ing PLUTO (Mignon e et al.ll2007h . a modular, Godunov-type code 
for astrophysical plasmas. The code provides a multiphysics, 
multi algorithm modular environment particularly oriented to- 
ward the treatment of astrophysical flows in the presence of 
discontinuities such as in the case treated here. The code was 
designed to make efficient use of massively parallel comput- 
ers using the message-passing interface (MPI) library for in- 
terprocessor communications. The MHD equations are solved 
using the MHD module available in pluto which is based on 
the Harten-Lax- van Le er Discontinuities (HLLD) approximate 
Riema nn solver jMiyoshi & Kusano 2005). Mivoshi & Kusano 
(120051) have shown that the HLLD algorithm can exactly solve 
isolated discontinuities formed in the MHD system and, there- 
fore, the adopted scheme is particularly appropriate to describe 
accretion shocks. The evolution of the magnetic field is carried 
out us ing the constrained transport method of Bal sara & Spiced 
d 19991) that maintains the solenoidal condition at machine accu- 
racy. PLUTO includes optically thin radiative losses in a fractional 
step formalism (Mignone et al. 2007), which preserve 2"'' time 
accuracy, being the advection and source steps at least 2""' order 
accurate; the radiative losses A values are computed at the tem- 
perature of interest using a table lookup/interpolation method. 
The thermal conduction is solved with an explicit scheme that 
adds the parabolic contributions to the upwind fluxes computed 
at cell interfaces (Mignone et al. 2007). Such a scheme is subject 
to a rather restrictive stability condition (i.e. At < (Ax)^/(2ri), 
where rj is the maximum diffusion coefficient), being the thermal 
conduction timescale generally shorter than the dynami cal one 
(e.g. Hujeirat & Camenzind 2000; ,Huieirat,2005; .Orlando et al.l 
I2005i 120081) . 

^ Note that the stream velocity is Ms„o = -580 km s"' when the stream 
impacts onto the chromosphere, due to gravity. 



The symmetry of the problem allows us to solve the MHD 
equations in half of the spatial domain with the stream axis co- 
incident with the z axis. The 2D cylindrical (r,z) mesh extends 
between and 1 .2 x 10'" cm in the r direction and between and 
2.4 X 10'° cm in the z direction and consists of a uniform grid 
with 512 X 1024 grid points. Additional runs were done with se- 
tups identical to those used for runs By-01, By-10, and By-50 
but with higher resolution (1024 x 2048 grid points; see Tab. 
[U to evaluate the effect of spatial resolution (see Sect. 13.4b ; the 
simulations with higher resolution cover a time interval of about 
1000 s. 

We use fixed boundary conditions at the lower (z = 0) bound- 
ary, imposing zero material and heat flux across the boundary. 
With this condition, matter may progressively accumulate at the 
base of the chromosphere; we have estimated that this effect may 
become significant on timescales a factor 100 longer than than 
those explored by our simulations. Axisymmetric boundary con- 
ditionflare imposed at r = (i.e. along the symmetry axis of the 
problem), a constant inflow at the upper boundary (z - 2.4x10'" 
cm) for r < rstr, and free outflo\\0 elsewhere. 

3. Results 

3. 1 . One-dimensional reference model 

The results of our ID reference hydrod ynamic simu l ation, HD- 
ID, are analogous to those discussed bv lSacco et al.l ( l2008h : the 
impact of the accretion stream onto the stellar chromosphere 
generates a reverse shock which propagates through the accre- 
tion column, producing a hot slab. According to Sacco et al., the 
expected temperature, Tsiab, and m aximum thickness, Dsiab, of 
the slab, in the strong shock limit dZel'Dovich & Raizeij|l96'7l) 
are: 

T,i,^^-^ul^«5MK, (10) 
iZ Kb 

Dsiab ~ 4.2 X 10^— t]'1 a; 2.5 X 10*^ cm , (11) 
nsti-o ''""' 

where is the Boltzmann constant. 

Fig. m shows the time-space plot of the temperature evolu- 
tion for run HD-ID. The base of the hot slab penetrates the 
chromosphere (the dashed line in Fig. |4] marks the initial po- 
sition of the transition region between the chromosphere and 
the corona) down to the position at which the ram pressure, 
^ram - Psti<)Ms,jO' °f post-sh ock plasma equals the thermal 
pressure of the chromosphere (iSacco et al.l 120081) . As evident 
from the figure, the shock front is not steady and the amplitude 
and period of the pulses rapidly reach stationary values as the 
initial transient disappears. The shock position oscillates with a 
period of « 600 s due to strong radiative cooling at the base of 
the slab, which robs the post-shock plasma of pressure support, 
causing the material above the c ooled layer to coll apse back be- 
fore the slab expands again (see ISacco et al.ll2008l for a detailed 
description of the system evolution). The maximum thickness of 
the slab is Dsiab ~ 2.5 x 10^ cm, in agreement with the predic- 
tion (see Eq.fTTTi. The post-shock plasma reaches the temperature 
T^siab ~ 5 - 7 MK during the expansion phase, and Tsiab ~ 2 - 4 
MK during the cooling phase. 



■* Variables are symmetrized across the boundary and normal com- 
ponents and angular components of vector fields («, B) flip signs. 
^ Set zero gradients across the boundary. 
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Fig. 4. Time-space plot of the temperature evolution for run HD- 
ID. The spatial extent of the shock lies in the vertical direction at 
any instant in time. The dashed line marks the initial position of 
the transition region between the chromosphere and the corona. 



As already mentioned, the evolution of the accretion shock 
described by the ID reference model is appropriate if the mag- 
netic field lines can be considered rigid to any force exerted by 
the accreting plasma flowing along the lines (yS «; 1). In case of 
plasma [3 values around 1 or even higher, the hot slab is expected 
to be only partially confined by the magnetic field and flow can 
occur sideways because of the strong pressure of the post-shock 
plasma and may ultimately perturb the dynamics of the shock 
itself. Moreover, 2D or 3D structures leading to strongly cool- 
ing zones are expected to develop in the post-shock plasma if 
the latte r is characterized by /? > 1. In this case. [Sutherland et al.l 
(I2003ah showed that the evolution of ID and 2D radiative shocks 
may be significantly different. 



3.2. Two-dimensional radiative MHD sinock model 

The cases of /3 > 1 are described by our 2D model. The main 
differences on the evolution from the ID model are expected to 
occur at the border of the stream where the shock-heated plasma 
might be not confined efficiently by the magnetic field; so we 
focus our attention there. 

Fig. |5] shows the spatial distribution of temperature (on the 
left) and plasma p (on the right) for runs By-01, By- 10, and By- 
50 at time t = 530 s (at early stage of evolution). Movies show- 
ing the complete evolution of 2D spatial distributions of mass 
density (on the left) and temperature (on the right), in log scale, 
for runs By-01, By- 10, and By-50 are provided as on-line ma- 
terial. The accretion flow follows the magnetic field lines and 
impacts onto the chromosphere, forming a hot slab at the base 
of the stream with temperatures x 5 MK and /3 > 1 . In all the 
cases, the slab is rooted down in the chromosphere, where the 
thermal pressure equals the ram pressure, and part (» 1 /3) of the 
shock column is buried under a column of optically thick ma- 
terial and may suffer significant absorption. In runs By-01 and 
By- 10, the dense hot plasma behind the shock front causes a 
pressure driven flow parallel to the stellar surface, expelling ac- 
creted material sideways (see upper and middle panels in Fig.|5j 
see also on-line movies for runs By-01 and By-10). This out- 
flow strongly perturbs the shock dynamics and is absent in the 
ID shock model. As expected, this feature determines the main 
differences between these 2D and the ID simulations. 



In run By-01, the magnetic field is too weak to confine the 
post-shock plasma (Ji ~ 10"^ at the border of the slab) and a 
conspicuous amount of material continuously escapes from the 
border of the accretion column at its base where the flow im- 
pacts onto the stellar surface. The maximum escape velocity is 
comparable to the free-fall velocity, MstrO, and this outflow acts 
as an additional cooling mechanism. The resulting outflow ad- 
vects and stretches the magnetic field lines (see upper panels in 
Fig. HI), taking the material away from the accretion column and 
strongly perturbing the stellar atmosphere even at several stream 
radii. As a result of the outflow, a strong component of B per- 
pendicular to the stream velocity (B, x 0.3|B()|) appears in the 
post- shock reg i on (se e upper panels in Fig.|5]). As discussed by 
Tot h & Draind (Il993h . the presence of even a small magnetic 
field perpendicular to the flow can stabilize the overstable os- 
cillations. In fact, at variance with the force due to gas pres- 
sure , the Lorentz force i s not affected by cooling processes (see 
also lHuieirat & Papaloizou 1998 ) and this mechanism may con- 
tribute to stabiUze the shock oscillations (see the on-line movie). 

In the intermediate run By-10, the magnetic field is trapped 
at the head of the escaped material, leading to a continuous in- 
crease of the magnetic pressure and field tension there. As a re- 
sult, the escaped material is kept close to the accretion column by 
the magnetic field (at variance with the By-01 case) and, even- 
tually, may perturb the stream itself. In fact, the escaped mate- 
rial accumulates around the accretion column, forming a grow- 
ing sheath of turbulent material gradually enveloping the stream. 
The magnetic pressure and field tension increase at the interface 
between the ejected material and the surrounding medium, push- 
ing on the material and forcing it to plunge into the stream after 
X 1.3 ks. Fig. |6] shows the temperature and mass density dis- 
tributions at time f = 1.8 ks, when the expelled material has 
already entered into the stream and has strongly perturbed the 
accretion column (see also the on-line movie to follow the com- 
plete evolution). As a result of the stream perturbation, the hot 
slab may temporarily disappear altogether as shown in Fig.|6l In 
this phase, the region of impact of the stream onto the chromo- 
sphere is characterized by a rather complex structure with knots 
and filaments of material and may involve possible mixing of 
plasma of the surrounding corona with accretion material. 

In the By-50 case, the post-shock plasma is confined effi- 
ciently by the magnetic field and no outflow of accreted mate- 
rial forms (see lower panels in Fig. [S]). In particular, /3 ~ 10""* 
and the magnetic pressure Pb ~ lOOf at the stream border, 
where Pram = pu^ is the ram pressure. The 2D shock, there- 
fore, evolves similarly to the ID overstable shock simulation, 
alternating phases of expansion and collapse of the post-shock 
region (see the on-line movie). The maximum thickness of the 
slab is Dsiab ~ 1.4 x 10^ cm, i.e. less than in the ID case by 
a factor =a 1.8. This re sult is analogous to that described by 
[Sutherland et al.l (l2003al) in the case of 2D hydrodynamic radia- 
tive shocks and is explained as due to the formation of denser 
knots of more rapidly cooling gas. Note that, at varia nce with our 
MHD s imulations, the hydrodynamic model of Sutherland et al.l 
(I2003al) does not include the thermal conduction. In our sim- 
ulations, the thermal conduction acts as an additional cooling 
mechanism of the hot slab, draining energy from the shock- 
heated plasma to the chromosphere, and partially contrasts the 
radiative cooling. Depending on the temperature of the shock 
and on the plasma /3, therefore, the thermal conduction may 
lead to significant differences between our results and those of 
Sutherland et al. (2003a). 

Fig. [T] shows snapshots of the evolution of temperature dis- 
tribution (in linear scale) in run By-50. The post-shock region 
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gets larger (f - 1.05 ks), reaches the maximum extension 
it - 1.15 ks), and collapses (f = 1.25 ks); the cycle then re- 
peats and the slab begins to get larger again (f = 1.35 ks). 
Although the magnetic field is strong enough to confine the post- 
shock plasma, the plasma fi is slightly larger than 1 inside the 
slab and the ID approxima tion is not valid there. A s in the 2D 
hydrodynamic simulations ( [Sutherland et al J [200331) . therefore, 
complex 2D cooling structures, including knots and filaments 
of dense material, form there due to the thermal instability of 
the post-shock plasma (see t - 1.15 ks in Fig. |7]l. As discussed 



by [Sutherland et al.l (l2003ab . these 2D complex structures lead 
to zones cooling more efficiently than those in ID models (by 
virtue of the increased cold-hot gas boundary) and, therefore, the 
amplitude of oscillations is expected to be reduced. It is worth 
noting that the similarity of our results with those of the hydro- 
dynamic models is due to the fact that, in run By-50, j6 > 1 in the 
slab and 2D cooling structures can form. In the case of yS «; 1 
everywhere in the slab, we expect that the stream can be consid- 
ered as a bundle of independent fibrils (each of them describable 
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in terms of ID models), and the 2D MHD simulations would 
produce the same results of ID models. 

We also study the global time evolution of 2D shocks by de- 
riving time-space plots of the temperature evolution analogous to 
that derived for the ID reference simulation (in order to compare 
directly our ID and 2D results). From the 2D spatial distributions 
of temperature and mass density, we first derive profiles of tem- 
perature along the z-axis, by averaging the emission-measure- 
weighted temperature along the r-axis for r < 4 x 10'' cm (i.e. 
80% of the stream radius); then, from these profiles, we derive 
the time-space plots of the average temperature evolution. 

Fig. [8] shows the results for runs By-01, By-10, and By-50, 
and can be compared with Fig. |4] In all the cases, the ampli- 
tude of the shock oscillations is smaller than that observed in the 
ID reference model. In runs By-01 and By-10, the evolution of 
the 2D shock is markedly different from that of the ID shock 
(compare Fig. |4] with top and middle panels in Fig. [SJ: in run 
By-01, the oscillations of the shock are stabilized after 2.5 ks 
and the solution approaches a quasi-stationary state; in run By- 
10, the oscillations appear chaotic without an evident periodicity 
(at least in the time lapse explored here). As discussed before, 
the stabilization of the shock oscillations in the former case is 
due t o a small magnetic fi eld component perpendicular to the 
flow dToth & Draine|[T993h . whereas the chaotic oscillations in 
the latter case are due to the continuous stream perturbation by 
the material ejected sideways at the stream base (see Fig|6]l. In 
the run By-50, the shock evolves somewhat similarly to the ID 
overstable shock simulation, showing several expansion and col- 
lapse complete cycles of the post-shock region (see lower panel 



in Fig.[8]l. However, as already discussed before, some important 
differences arises: by comparing By-50 with HD-ID, the oscil- 
lations are reduced in amplitude by a factor ~ 1.8 and occur at 
higher frequency (period Psh ~ 300 s). 

3.3. Effect of stream parameters 

The details of the shock evolution described in this paper depend 
on the model parameters adopted. In particular, the temperature 
of the shock-heated plasma is determined by the free-fall ve- 
locity with which the plasma impacts onto the star (EqfTOli; the 
stand-off height of the hot slab generated by the impact depends 
on the velocity and density of the stream (see Eq.fTTI). Assuming 
a typical value for the free-fall velocity of 400 - 500 km s"' 
(leading to temperatures of » 3 - 5 MK, as deduced from obser- 
vations), the maximum thickness of the slab is determined only 
by the density of the stream: the heavier the stream, the thinner 
is expected the slab. Also the sinking of the stream down in the 
chromosphere depend on the stream density and velocity at im- 
pact. In fact, we find that the slab penetrates the chromosphere 
down to the position at which the ram pressure, PstroMstiO' '-^^ 
post-shock plas ma equals the ther mal pressure of the chromo- 
sphere (see also Sacco et al.ll2008 ). As a result, heavier streams 
are expected to sink more deeply into the chromosphere. 

In the case of /3 > 1, the shock evolution is expected to be 
influenced also by the stream radius, rjtr. In fact, as shown by our 
simulations, the complex plasma dynamics close to the border of 
the stream (e.g. the generation of outflows of accreted plasma; 
see runs By-01 and By-10) may affect the shock evolution in the 
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inner portion of the slab. The shock evolution is expected to be 
modified in the whole slab if the oscillation period of the shock, 
Psh, is larger than the dynamical response time of the post-shock 
region, which can be approximated as the sound crossing time of 
half slab: Tdyn ~ rsu/cs, where Cj is the isothermal sound speed. 
In the cases discussed in this paper, the oscillation period is Psh ~ 
300 s, and the sound crossing time of the slab is Tdyn ~ 200 s; as 
shown by our simulations, the shock evolution is modified in the 
whole slab. Assuming an accretion stream with twice the radius 
considered here, we find Tdyn ~ 400 s, and the region at the 
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Fig. 8. As in Fig.|4]for the 2D simulations By-01 (top), By-10 
(middle) and By-50 (bottom). The blue curve in the bottom panel 
marks the shock position derived in the 1 D reference model, HD- 
ID. 



center of the slab should not be affected by the plasma dynamics 
at the stream border. In this case, we expect quasi-periodic shock 
oscillations (analogous to those described by run By-50) at the 
center of the stream and a strongly perturbed shock at the stream 
border. 



3.4. Effect of spatial resolution 

In problems involving radiative cooling, the spatial resolution 
and the numerical diffusion play an important role in determin- 
ing the accuracy with which the dynamics of the system is de- 
scribed. In particular, in the case of radiative shocks, we expect 
that the details of the plasma radiative cooling depend on the 
numerical resolution: a higher resolution may lead to different 
peak density and hence influence the cooling efficiency of the 
gas, and therefore the amplitude and frequency of shock oscilla- 
tions. In the simulations presented here, the thermal conduction 
partially contrasts the radiative cooling all eviating the proble m 
of numerical resolution (see, for instance, Orlando et al.ll2008h . 

In order to check if our adopted resolution is sufficient to 
capture the basic shock evolution over the time interval consid- 
ered, we repeated simulations By-01, By-10, and By-50 but with 
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Fig. 9. Evolution of the averaged shock-front location for runs 
By-50 (solid blue Hne) and By-50-HR (dashed red line). The 
dotted hne marks the initial position of the transition region be- 
tween the chromosphere and the corona. 

twice the spatial resolutior0 (runs By-Ol-HR, By- 10-HR, and 
By-50-HR). The efficiency of radiative cooling is expected to be 
the largest in the model with |Bol = 50 G, because there is no loss 
of material through sideways outflows. In fact, quasi-periodic os- 
cillations of the accretion shock occur (see Fig.lSJ and make this 
case adequate for a comparison of different spatial resolutions. 
Since this case is one of the most demanding for resolution, it 
can be considered a worst case comparison of convergence. 

Fig. |9] compares the evolution of the averaged shock-front 
location for runs By-50 and By-50-HR. The main differences 
between the simulations appear in the first bounce which starts 
earlier and is larger by a factor 2 in By-50-HR than in By-50. 
The first expansion of the hot slab occurs after that the stream 
has penetrated the chromosphere down to the position where the 
ram pressure of the shock-heated material equals the thermal 
pressure of the chromosphere. The first bounce is a transient, 
therefore, and in fact its amplitude and width are different from 
those of subsequent bounces in all the simulations examined. 
On the other hand, Fig.|9]clearly shows that, apart from the first 
bounce, in general, the results of the two simulations agree quite 
well, showing differences < 10%. We expect, therefore, that in- 
creasing the spatial resolution, the results of our simulations may 
slightly change quantitatively (e.g. the amplitude of oscillations) 
but not qualitatively. 

4. Discussion 

4.1. Distribution of emission measure vs. temperature of the 
sliock-lieated plasma 

As discussed in the introduction, there is a growing consen- 
sus that two distinct plasma components contribute to the X-ray 
emission of CTTSs: the stellar corona and the accretion shocks. 
This idea has been challenged recently by Argiroffi et al.l (l2009l) . 
who compared the distributions of emission measure, EM(r), 
of two CTTSs with evidence of X-ray emitting dense plasma 
(MP Mus and TW Hya) with that of a star (TWA 5) with no ev- 
idence of accretion (i.e. only the coronal component is present) 

* Note that, in the case of simulations with higher spatial resolution 
(very CPU time consuming), the time interval covered was ~ 1 ks, in- 
stead of 3 ks, in order to reduce the computational cost. 



and with the EM(r) derived from a ID hydrodynamic model of 
accretion s hocks (i .e. only the shock-heated plasma component 
is present; ISacco e t al. 2008). They proved that the EM(T) of 
MP Mus and TW Hya can be naturally interpreted as due to a 
coronal component, dominating at temperatures T > 5 MK, plus 
a shock-heated plasma component, dominating at T < 5 MK. 

In case of shocks with yS > 1, our simulations show that the 
distributions of temperature and density at the base of the ac- 
cretion column can be rather complex (see, for instance, Fig.|7j. 
It is interesting, therefore, to investigate how the EM(r) of the 
shock-heated plasma changes with /3 and if it is possible to derive 
a diagnostics of the plasma-/? in the post-shock region. 

In order to derive the EM(r) distribution of the accretion 
region from the models, we first recover the 3D spatial distribu- 
tions of density and temperature by rotating the coiTesponding 
2D distributions around the symmetry z axis (r - 0). The emis- 
sion measure in the jth domain cell is calculated as emj - Vj, 

where is the particle number density in the cell, and Vj is the 
cell volume. The EM(r) distribution is then derived by binning 
the emission measure values into slots of temperature; the range 
of temperature [5 < log r(K) < 8] is divided into 15 bins, all 
equal on a logarithmic scale (A log T - 0.2). 

Fig.[TO]shows the EM(r) distributions averaged over 3 ks for 
runs By-01, By-10, and By-50, together with the average EM(r) 
derived from our ID reference model HD-ID (blue dashed 
lines). The figure shows the EM(r) distributions of plasma with 
density tin > 10'' cm"^ (black) and 10'" < ne < 10" cm"^ 
(red). The 2D MHD models have been normalized in order to 
have an X-ray luminosity Lx ~ 4x10"'' erg in the band [0.5-8.0] 
keV, in agreement with the luminosity derived from the low- 
temperatu re (log T < 6.7) port ion of the EM(r) distribution of 
MP Mus (lArgiroffi et al.ll2009h . This is obtained assuming that 
w 10 accretion streams similar to that modeled here are present 
simultaneously. The ID model has been normalized in order to 
match the EM peak in By-50. Inspecting Fig.[TO]we note that: i) 
in all the cases, the EM(r) has a peak at ~ 5 MK and a shape 
compatible with those observed in MP Mus and TW H ya and 
attributed to shock-heated material (e.g. lArgiroffi et al.ll2009l) ; 
ii) most of the X-ray emission arises from shock-heated plasma 
with density > 10" cm"^ regardless of the /3 (i.e. the most dense 
component of the post-shock region dominates); iii) the slope of 
the ascending branch of the EM(r) distribution is comparable in 
runs By-50 and HD-ID, and gets steeper for decreasing values 
of /3. The time-averaged X-ray luminositjQ derived from these 
EM(r) distributions ranges between 5 x 10^^ erg (run By-10) 
and 8 x 10^^ erg (run By-01), having the X-ray luminosity only 
a weak dependence on the plasma /3. 

Our model supports, therefore, the idea that dense shock- 
heated plasma may contribute significantly to the low tempera- 
ture portion of the EM(r) distributions of CTTSs regardless of 
the /? value. We also suggest that the shape of the EM(r) might 
be used as a diagnostics of /3 in the post-shock region, if the coro- 
nal contribution to the low temperature tail can be neglected. 

4.2. Mass accretion rates 

Time-dependent models of radiative accretion shocks provide a 
strong theoretical support to the hypothesis that soft X-ray emis- 
sion from CTTSs arises from shocks due to the impact of the ac- 

' The synthetic X-ray spectra in the band [0. 5 - 8.0] keV have been 
derived with the PINTofALE spectral code ( Kashvap & Drake! [2000h 
with the APED VI. 3 atomic line database ( Smit h et alj200lh . assuming 
the metal abundances adopted in the whole paper, namely 0.5. 



10 



S. Orlando et al.: X-ray emitting MHD accretion shocks in CTTSs 



53 



§ 52 



LU 

o 



51 r. 



50 

49 
53 



§ 52 



LU 
O 



LU 
O 



51 

50 

49 
53 

52 

51 

50 
49 



10^° < Hh < 10^^ cnn-3 
nn > 10^^ cm"^ 
HD-1D 



By-01 



By-10 



By-50 



5.5 



6.0 



6.5 
Log T [K] 



7.0 



7.5 



Fig. 10. Distributions of emission measure vs. temperature, av- 
eraged over 3 ks, for runs By-01, By-10, and By-50. Black (red) 
lines mark the average EM(r) distributions of plasma with den- 
sity riu > 10" cm"-' (10'" < hh < 10" cm"-'); dashed blue lines 
mark the average EM(r) distribution for the ID reference model 
HD-ID. 



cretion columns onto the stellar surface. However, several points 
remain still unclear. Among these, the most puzzling is probably 
the fact that the mass accretion rates derived from X-rays, Mx, 
are consistently lower by one or more orders of magnitude than 
the correspondi ng M v alues de rived from UV /optical/ NlR ob- 
serva t ions (e.g.lDrakdb OOS: Sc hmitt et al]|200 5; Gunth er et all 
l2007t lArgiroffi et aiTl2009i CuiTan et al. 2009, in preparation). 
We have here the opportunity to discuss the problem of the ac- 
cretion rate, in particular focusing on the cases with yS > 1 in the 
post-shock region. 

The model parameters adopted in this paper describe a 
stream with an accretion rate of ^ 9 x lO^'^M© yr"'. According 
with the discussion in Sect. 14.11 ^ 10 streams are needed to 
match the soft X-ray luminosity of MP Mus. In this case the 
accretion rate is ^ 9 x 1O""M0 yr"', which is in nice agree- 
ment with that deduced from observ ations, namely Mx = 5 - 8 x 



10~"Mo yr"' dArgiroffi et alj|2007h . Alternatively, the observed 
Mx may be reproduced by our model if the accretion stream has 



a larger cross section (with radius r^a- ~ 10'" cm). In this case, 
we expect some changes to the dynamics of the shock-heated 
plasma as described in Sect. 13.31 

On the other hand, the mass accretion rate of MP Mus, as 
deduced from optica l observations, is Mopi ~ 3 x lO^^M© yr"' 
dArgiroffi et al.l2009l) . exceeding by more than one order of mag- 
nitude the value obtained from X-rays. Analogous discrepancies 
are found in all CTTSs for which it is possible to derive Mx 
(see also Curran et al. 2009, in preparation). The discrepancy 
might be reconciled if Mx values are underestimated due, for 
instance, to absorption from optically thick plasma. However, 
as we explain below, even assuming that the absorption can ac- 
count for the observed M discrepancy, the idea that the same 
streams determine both Mopt and Mx has to be discarded. In 
fact, assuming the accretion parameters adopted here, we derive 
that a! 300 streams must be present to match Mopt or the stream 
should have a cross section a 2 x 10^^ cm^, implying that, in 
both cases, ^ 20 % of the stellar surface should be involved in 
accretion. None of these hypothesis is realistic, being the sur- 
face filling facto r of hot spo ts resulting from accretion up to a 
few percent (e.g.lBouviei^& Bertout 1989; Hartmann & K en^onl 
1990; Bouvieretal. 1993; Gullbring 1994; IKenvon etalj|1994l; 
iDrake et al.ii2009i and references therein). 

The above inconsistency can be removed if the mass density 
of accretion streams is larger than assumed here; for instance, for 
a density of the stream «stro = 5 x 10'^ cm"-* (a factor 50 larger 
than modeled), the mass accretion rate is a; 5 x 10" '"Mq yr"' and 
few streams are needed to match Mopt- In this case, however, the 
model cannot explain the lower value s of density derived from 
X-ray observations (nstiO ~ 10" cm"^ ; lArgirofR et ani2007h . 

A possible solution to remove the discrepancy in the M val- 
ues (and in the WstrO values) is that few accretion streams charac- 
terized by different mass density values are present at the same 
time: those with density x 10" cm"^ would produce shocks that 
are visible in the X-ray band, leading to the observed Mx; those 
with higher densities would produce shocks not visible in the X- 
ray band and leading to the observed Mopt (being the dominant 
component in the UV/optical/NlR bands). The reason of the dif- 
ferent visibility could be due to local absorption of the X-ray 
emission, as explained in the following. 

Our simulations show that the accretion stream penetrates 
the chromosphere down to the position at which the ram pressure 
of the post-shock plasma equals the thermal pressure of the chro- 
mosphere (see also Sacco et al. 2008). Part of the shock-heated 
plasma is buried down in the chromosphere and is expected to be 
obscured by significant absorption from optically thick plasma 
(see also Sacco et al. 2009, in preparation). In the simulations 
presented here, this portion is 1/3 of the hot slab (/isink ~ 5 x 10^ 
cm; see Fig. [8j, and most of the post-shock plasma (above the 
chromosphere) is expected to be visible with minimum absorp- 
tion. Instead, in denser streams, the shock column is buried more 
deeply in the chromosphere, due to the larger ram pressure, and 
its maximum thickness is smaller, according to Eq.[TT](see Sect. 
13.31 ). Assuming the accretion parameters adopted in this paper 
but with Mstio = 5 X 10'^ cm"-', Eq.[TT] gives Aiab 5 x 10^ 
cm, that is much smaller than the expected sinking of the stream 
in the chromosphere (/isink ~ 7 x 10*^ cm). In this case, the 
post-shock column is buried under a hydrogen column density 
Nh = nuhsink ~ 2 X 10^^ cm"^ and the photoelectric absorption 
of the Ovii triplet (i.e. the lines commonly used to trace the ac- 
cretion in the X-ray band) is given by exp[-croviiA^H] ~ 10"^. 
We conclude, therefore, that heavy streams may produce X-ray 
emitting shocks that, however, are hardly visible in X-rays, being 
buried too deeply in the chromosphere. 
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4.3. Variability of X-ray emission from sfiock-lieated plasma 

Another important point in the study of accretion shocks in 
CTTSs is the periodic variability of X-ray emission, due to 
the quasi-periodic oscillations of the shock position induced by 
cooling, predicted by time-dependent ID models. However, in 
the only case analyzed up to date, namely TW Hya, no ev- 
idence of periodic variations of soft X-ray emission (thought 
to arise p redominantly in an accretion shock) has been found 
jPrake et al. 2009). This result apparently contradicts the predic- 
tion of current ID models and Drake et al. suggested that these 
models might be too simple to explain the 3D shock structure. 

On the other hand, quasi-periodic shock oscillations are ex- 
pected if the accretion stream is homogeneous and constant (no 
variations of mass density and velocity). However, there is sub- 
stantial observational evid ence that the s t reams are clumped 
and inhomogeneous (e^ iGuUbring et all [19961: ISafieJ [19981 
iBouvier et al. 2003, 2003). In these conditions, periodic shock 
oscillations are expected to be hardly observable. 

The 2D MHD simulations presented here show that, even as- 
suming constant stream parameters, periodic oscillations are not 
expected if/? » 1 in the post-shock region (runs By-01 and By- 
10). In these cases, in fact, the time-space plots of temperature 
evolution (top and middle panels in Fig.[8]l predict that the oscil- 
lations may be rapidly dumped, approaching a quasi-stationary 
state with no significant variations of the shock position, or the 
variability may be chaotic (with no obvious periodicity) due to 
strong perturbation of the stream by the accreted material ejected 
sideways. In none of these cases, therefore, we expect to observe 
periodic modulation in the X-ray emission. 

On the other extreme, for shocks with /3 <*c 1, we expect 
that the single accretion stream is structured in several fibrils, 
each independent from the others due to the strong magnetic 
field which prevents mass and energy exchange across magnetic 
field lines. Time-dependent ID models describe one of these fib- 
rils. Being independent from each other, the fibrils can be char- 
acterized by slightly different mass density and velocity, which 
would result in different instability periods, a s well as by ran- 
dom phases of the oscillations. Dra ke et akl (|2009l). assuming 
that ID models describe a single stream, derived the number of 
streams needed to account for the absence of periodic variabil- 
ity in TW Hya and concluded that this number contrasts with 
the presence of conspicuous rotatio nally m odulated surface flux 
with small filling factor Following iDrake et al. (2009) but con- 
sidering the fibrils instead of the streams, we may argue that an 
accretion stream consisting of 200 - 300 different fibrils with 
radius rfjbr ~ 10^ cm and with different instability periods and 
random phases would produce a signal pulsed at a level of less 
than 5% as measured in TW Hya. 

The intermediate situation is for shocks with around 1 . In 
this case, our run By-50 shows quasi-periodic oscillations of the 
shock position (see bottom panel in Fig.|8]l and predicts periodic 
variations of the X-ray emission arising from a single stream. In 
fact, in this case, the magnetic field is strong enough to confine 
the shock-heated plasma, but it is too weak to consider valid the 
ID approximation inside the slab that cannot be described as a 
bundle of fibrils. As a result, the shock o scillates coherently in 
the slab. As noted by Drake et al.l (|2009), in this case it is not 
possible to reproduce the absence of periodic modulation of X- 
ray emission observed in TW Hya, and we conclude that shocks 
with jS 1 do not occur in this star. 

It is worth emphasizing that further investigation is needed 
to understand how much the intermediate case described by run 
By-50 is frequent and observable: an intensive simulation cam- 



paign is needed to assess the range of /3 values leading to streams 
with quasi-periodic oscillations: also a systematic analysis of the 
variability of soft X-ray emission should be performed, consid- 
ering a complete sample of CTTSs with evidence of X-ray emit- 
ting accretion shocks, to assess if and when periodic variability 
is observed. 

4.4. Effects of accretion shocks on the surrounding stellar 
atmosphere 

In the case of shocks with y6 » 1, our simulations predict that the 
stellar atmosphere around the region of impact of the stream can 
be strongly perturbed by the impact, leading to the generation of 
MHD waves and plasma motion parallel to the stellar surface: 
the larger is the plasma /3 in the post-shock region, the larger is 
the perturbation of the atmosphere around the shocked slab. The 
resulting ejected flow also advects the weak magnetic field such 
that the conditions for ideal MHD may break down, magnetic 
reconnection may be possible and eventually releases the stored 
energy from the magnetic field (not described however by our 
model that does not include resistivity effects). 

A possible effect of the perturbation of the stellar atmosphere 
is that shocks wi t h P » 1 may contribute to the stellar outflow. 
In fact, ICranme 3 (l2008h suggested that the MHD waves and the 
material ejected from the stream (as in our runs By-01 and By- 
10) may trigger stellar outflow and proposed a theoretical model 
of accretion-driven winds in CTTSs (see, also, Cranmer 200^. 
His model originates from a description of the coronal heating 
and wind acceleration in the Sun and includes a source of wave 
energy driven by the impact of accretion streams onto the stellar 
surface (in addition to the convection-driven MHD turbulence 
which dominates in the solar case). The author found that this 
added energy seems to be enough to produce T Tauri-like mass 
loss rates. It would be interesting to assess how the different 
plasma-/? cases discussed in this paper contribute to the added 
energy. 

4.5. Limits of the model 

Our simulations were carried out in 2D cylindrical geometry, 
implying that all quantities are cyclic on the coordinate (p. This 
choice is expected to affect some details of the simulations but 
not our main conclusions. In particular, adopting a 3D Cartesian 
geometry, the simulations would provide an additional degree of 
freedom for hydrodynamic and thermal instabilities, increasing 
the complexity of cooling zones and, therefore, the cooling ef- 
ficiency of the plasma in the post-shock region. As a result, in 
3D simulations, the amplitude of the shock oscillations might 
be slightly smaller and the frequency slightly higher than that 
observed in 2D simulations. 

As discussed in Sect. 13.31 some details of our simulations 
depend on the choice of the model parameters. For instance, the 
temperature of shock-heated plasma, the stand-off height of the 
hot slab as well as the sinking of the stream down in the chromo- 
sphere depend on the stream density and velocity at impact. The 
accretion parameters adopted here originate from the values de- 
rived bv [Argiroffi et al.l (120071) for MP Mus; the cases presented 
here are representative of a regime in which the shock-heated 
plasma has a temperature x 5 MK and part of the slab is above 
the chromosphere, making it observable with minimum absorp- 
tion. The range of magnetic field strength considered has been 
chosen in order to explore shocks with /? > 1 . Nevertheless, our 
results undoubtedly show that the stability and dynamics of ac- 
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cretion shocks strongly depends on /? and that a variety of phe- 
nomena, not described by ID models, arises, including the gen- 
eration of plasma motion parallel to the stellar surface and MHD 
waves. 

It is worth noting that our model does not include magnetic 
resistivity effects. The complex shock evolution described by our 
simulations show that violent plasma flows may advect the mag- 
netic field such that the conditions for ideal MHD may break 
down and magnetic reconnection may possibly occur (not de- 
scribed by our simulations). In particular, field lines at the stream 
border are squeezed close together and may reconnect, result- 
ing in a change of the magnetic topology as well as in a release 
of magnetic energy. Since most of these events are expected to 
occur at the stream border, magnetic reconnection may play an 
important role in the dynamics and energetic of the ejected ac- 
creted plasma. For instance, the magnetic reconnection may trig- 
ger flares at the stream border outside the accretion column. This 
scenario is not in contradiction with Re ale (.2009.) , who showed 
that it is unlikely that flares can be triggered in an accreting flux 
tube. This issue deserves further investigation in future studies. 

5. Conclusion 

We investigated the stability and dynamics of accretion shocks 
in CTTSs, considering the case of y6 > 1 in the post-shock re- 
gion, through numerical MHD simulations. To our knowledge, 
the simulations presented here represent the first attempt to 
model 2D accretion shocks that simultaneously includes mag- 
netic fields, radiative cooling, and magnetic-field-oriented ther- 
mal conduction. Our findings lead to several conclusions: 

1 . In all the cases, a hot slab of shock-heated material is gener- 
ated at the base of the accretion column due to the impact of 
the stream with the chromosphere. In the case of shocks with 
y6 > 10, violent outflows of shock-heated material and, possi- 
bly of MHD waves, are generated at the border of the hot slab 
and they may perturb the surrounding stellar atmosphere. For 
shocks with (3 ~ I, the shock-heated plasma is efliciently 
confined by the magnetic field and no outflow forms. 

2. If the magnetic field is too weak to confine the shock-heated 
plasma but is strong enough to keep it close to the stream, the 
escaped accreted material may strongly perturb the accretion 
column, modifying the dynamics and stability of the shock 
itself. 

3. The accretion shocks are susceptible to radiative shock in- 
stability. The resulting shock oscillations strongly depend 
on the plasma for fi > 10, the oscillations may be 
rapidly dumped by the magnetic field, approaching a quasi- 
stationary state, or may be chaotic with no obvious period- 
icity due to perturbation of the stream induced by the post- 
shock plasma itself; for p around 1 the oscillations are quasi- 
periodic with amplitude smaller and frequency higher than 
those predicted by ID models. 

We also suggest that the M discrepancy observed in CTTSs 
may be solved if few streams with significantly different mass 
density are present: heavy streams are not visible in X-rays due 
to strong absorption and determine the M value deduced from 
observations in UV/optical/NIR bands; light streams can be ob- 
served in X-rays and determine the M value deduced from X- 
ray observations. As for the absence of periodic X-ray modula- 
tion du e to shock oscillations found in TW Hya by Drake et al. 
(I2009I) . we find that no periodic modulation of X-rays is ex- 
pected in cases of shocks with either y6 « 1 or/? » 1, whereas 



periodic modulation is found for fi around 1 , indeed a rather lim- 
ited set of possible cases. We interpret the absence of periodic 
X-ray modulation in TW Hya as evidence of the fact that, in this 
star, fi differs significantly from 1 in the accretion shocks or the 
streams are inhomogeneous. 
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